Compiled under R version 3.6.3 (2020-02-29)

WARNING: edit the working directory to your preferred folder.

This document details all analyses performed in R for the study:
Legendre, L. J., D. Rubilar-Rogers, G. M. Musser, S. N. Davis, R. A. Otero, A. O. Vargas, and J. A. Clarke. 2020. A giant soft-shelled egg from the Late Cretaceous of Antarctica. Nature 583, 411–414.

For more information regarding the study, datasets, and analyses, please refer to the Supplementary Information of this paper. If you have any additional questions, feel free to email me at .

Loading packages

library(AICcmodavg)
library(ape)
library(caper)
library(dplyr)
library(evobiR)
library(ggplot2)
library(MPSEM)
library(nlme)
library(nortest)
library(phylopath)
library(phytools)
library(RColorBrewer)

1 – Lepidosaur dataset

This part requires the use of dataset 1 (“Dataset1-lepidosaurs.txt”; see also Supplementary Table 1 in Legendre et al., 2020), and lepidosaur phylogenetic tree (“Lepidosaurtree.trees.nex”).

  • Loading the tree:
tree<-read.nexus("Lepidosaurtree.trees.nex")

Phylogenetic generalized least squares (PGLS)

Egg volume ~ SVL (snout-vent length)

  • Loading the dataset and removing species with missing data
data2<-read.table("Dataset1-lepidosaurs.txt", header=T); data2V<-subset(data2, !is.na(V)&!is.na(SVL))
pruned.tree<-drop.tip(tree, setdiff(tree$tip.label, data2V$Species))
data<-read.table("Dataset1-lepidosaurs.txt", header=T, row.names = "Species"); dataV<-subset(data, !is.na(V)&!is.na(SVL))
  • Best fit for alpha parameter in OU (Ornstein-Uhlenbeck) model
alpha <- seq(0, 1, 0.1)
fit <- list()
form<-log(V)~log(SVL)
for (i in seq_along(alpha)) {
  cor <- corMartins(alpha[i], phy = pruned.tree, fixed = T)
  fit[[i]] <- gls(form, correlation = cor, data = data, na.action=na.exclude, method = "ML")
}
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
sapply(fit, logLik)
##  [1]       Inf -325.1206 -325.9100 -326.0731 -326.1266 -326.1503 -326.1665
##  [8] -326.1806 -326.1939 -326.2064 -326.2182
  • Best fit for g parameter in EB (Early Burst) model
g <- seq(0.1, 1, 0.1)
fit <- list()
form <- log(V)~log(SVL)
for (i in seq_along(g)) {
  cor <- corBlomberg(g[i], phy = pruned.tree, fixed = T)
  fit[[i]] <- gls(form, correlation = cor, data = data, na.action=na.exclude, method = "ML")
}
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
sapply(fit, logLik)
##  [1] -322.9247 -331.3311 -337.6858 -341.8988 -344.8120 -346.9305 -348.5380
##  [8] -349.8001 -350.8185 -351.6589
  • Building PGLS models and running model selection using AICc
BM<-gls(log(V)~log(SVL), data=dataV, correlation=corBrownian(phy=pruned.tree), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OU<-gls(log(V)~log(SVL), data=dataV, correlation=corMartins(0.1, phy=pruned.tree, fixed=T), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
Lambda<-gls(log(V)~log(SVL), data=dataV, correlation=corPagel(1, phy=pruned.tree), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
EB<-gls(log(V)~log(SVL), data=dataV, correlation=corBlomberg(0.1, phy=pruned.tree, fixed=T), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OLS<-gls(log(V)~log(SVL), data=dataV, method="ML")

Cand.models = list()
Cand.models[[1]] = BM
Cand.models[[2]] = OU
Cand.models[[3]] = Lambda
Cand.models[[4]] = EB
Cand.models[[5]] = OLS

Modnames = paste(c("BM", "OU", "Lambda", "EB", "OLS"), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = T)
summary(Lambda)
## Generalized least squares fit by maximum likelihood
##   Model: log(V) ~ log(SVL) 
##   Data: dataV 
##        AIC      BIC    logLik
##   641.6771 655.6163 -316.8385
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.5664802 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 1.978195 0.5517311  3.585434   4e-04
## log(SVL)    1.114992 0.0559233 19.937869   0e+00
## 
##  Correlation: 
##          (Intr)
## log(SVL) -0.51 
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -3.0709107 -0.7544921 -0.1511364  0.4256780  2.0898411 
## 
## Residual standard error: 1.157687 
## Degrees of freedom: 241 total; 239 residual
  • Testing for normality
plot(Lambda, resid(., type="n")~fitted(.), main="Normalized Residuals v Fitted Values",abline=c(0,0))

res<-resid(Lambda, type="n")
par(mar=c(5.1,4.1,4.1,2.1))
qqnorm(res)
qqline(res)

lillie.test(residuals(Lambda))
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuals(Lambda)
## D = 0.053479, p-value = 0.09167
lillie.test(chol(solve(vcv(pruned.tree)))%*%residuals(Lambda))
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  chol(solve(vcv(pruned.tree))) %*% residuals(Lambda)
## D = 0.078349, p-value = 0.001089
  • Removing outliers and building a new model without them
which(res>4); which(res<(-2.5))
## named integer(0)
##      Chirindia_ewerbecki         Lycodon_striatus      Micrurus_mipartitus 
##                       89                      196                      222 
##         Ogmodon_vitianus Leptotyphlops_scutifrons    Myriopholis_blanfordi 
##                      224                      227                      228 
##    Indotyphlops_braminus 
##                      235
newdata<-dataV[-c(89,168,222),]
newdata2<-data2V[-c(89,168,222),]
newtree<-drop.tip(pruned.tree, setdiff(pruned.tree$tip.label, newdata2$Species))

Newlambda<-gls(log(V)~log(SVL), data=newdata, correlation=corPagel(1, phy=newtree), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
summary(Newlambda)
## Generalized least squares fit by maximum likelihood
##   Model: log(V) ~ log(SVL) 
##   Data: newdata 
##        AIC      BIC    logLik
##   615.7252 629.6143 -303.8626
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.3749586 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 1.773470 0.4386847  4.042699   1e-04
## log(SVL)    1.137279 0.0542297 20.971523   0e+00
## 
##  Correlation: 
##          (Intr)
## log(SVL) -0.618
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -3.4479913 -0.7458697 -0.0785325  0.6087246  2.5045329 
## 
## Residual standard error: 0.9939715 
## Degrees of freedom: 238 total; 236 residual
plot(Newlambda, resid(., type="n")~fitted(.), main="Normalized Residuals v Fitted Values",abline=c(0,0))

res<-resid(Newlambda, type="n")
par(mar=c(5.1,4.1,4.1,2.1))
qqnorm(res)
qqline(res)

lillie.test(residuals(Newlambda))
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuals(Newlambda)
## D = 0.057905, p-value = 0.05125
lillie.test(chol(solve(vcv(newtree)))%*%residuals(Newlambda))
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  chol(solve(vcv(newtree))) %*% residuals(Newlambda)
## D = 0.091664, p-value = 4.941e-05
  • Estimating the pseudo R-squared using caper
#newlambda:
datacomp<-comparative.data(phy=newtree, data=newdata2, names.col="Species")
pgls<-pgls(log(V)~log(SVL), data=datacomp, lambda="ML")
summary(pgls)
## 
## Call:
## pgls(formula = log(V) ~ log(SVL), data = datacomp, lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.107910 -0.039203 -0.005097  0.034535  0.147102 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.718
##    lower bound : 0.000, p = 1.8499e-05
##    upper bound : 1.000, p = 4.5223e-05
##    95.0% CI   : (0.413, 0.908)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) -0.41450    0.91690 -0.4521    0.6528    
## log(SVL)     1.63215    0.15212 10.7293 6.661e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06154 on 64 degrees of freedom
## Multiple R-squared: 0.6427,  Adjusted R-squared: 0.6371 
## F-statistic: 115.1 on 1 and 64 DF,  p-value: 6.121e-16
#lambda:
datacompl<-comparative.data(phy=pruned.tree, data=data2V, names.col="Species")
pgls<-pgls(log(V)~log(SVL), data=datacompl, lambda="ML")
summary(pgls)
## 
## Call:
## pgls(formula = log(V) ~ log(SVL), data = datacompl, lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.140589 -0.048249 -0.004924  0.031043  0.155446 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.715
##    lower bound : 0.000, p = 7.7636e-05
##    upper bound : 1.000, p = 0.00015352
##    95.0% CI   : (0.381, 0.918)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  0.19647    0.96534  0.2035    0.8393    
## log(SVL)     1.53068    0.16093  9.5115 5.418e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06683 on 66 degrees of freedom
## Multiple R-squared: 0.5782,  Adjusted R-squared: 0.5718 
## F-statistic: 90.47 on 1 and 66 DF,  p-value: 5.417e-14

Not much difference between the two models, assumptions of normality are not met in either of them – we use Lambda instead of Newlambda for the plot.

  • Plot of the regression
ggplot(dataV, aes(log(SVL), log(V), color=Suborder)) +
geom_point(size=4) +
  xlab("ln snout-vent length (mm)") +
  ylab("ln egg volume (mm3)") +
  geom_abline(intercept=Lambda$coefficients[1], slope=Lambda$coefficients[2], colour="lightblue", size=1.3) +
  theme(panel.background = element_rect(fill="black")) +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15),
        legend.text=element_text(size=15),
        legend.title=element_text(size=17)) +
        scale_colour_brewer("Suborder", palette="Accent")

  • Estimating the SVL of Antarcticoolithus using the allometric equation of PGLS, including it in the dataset, then in the plot
exp((-0.7990941/1.6708400)+(1/1.6708400)*log(5471405.789))
## [1] 6684.303
data[1,7]<-6684.303
dataV2<-subset(data, !is.na(V))

ggplot(dataV2, aes(log(SVL), log(V), color=Suborder)) +
  geom_point(size=4) +
  xlab("ln snout-vent length (mm)") +
  ylab("ln egg volume (mm3)") +
  geom_abline(intercept=Lambda$coefficients[1], slope=Lambda$coefficients[2], colour="lightblue") +
  theme(panel.background = element_rect(fill="black")) +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15),
        legend.text=element_text(size=15),
        legend.title=element_text(size=17)) +
  scale_colour_brewer("Suborder", palette="Accent")

Eggshell thickness ~ body mass is not detailed here, since the code is virtually identical – only variable names and regression coefficients for the plot and estimation are different.

  • Estimating the body mass of Antarcticoolithus using the allometric equation of PGLS, including it in the dataset, then in the plot
LambdaBM<-gls(log(V)~log(BM), data=dataV, correlation=corPagel(1, phy=pruned.tree), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
summary(LambdaBM)
## Generalized least squares fit by maximum likelihood
##   Model: log(V) ~ log(BM) 
##   Data: dataV 
##      AIC      BIC  logLik
##   524.34 538.2792 -258.17
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.3112388 
## 
## Coefficients:
##                Value  Std.Error  t-value p-value
## (Intercept) 5.260359 0.26352407 19.96159       0
## log(BM)     0.645885 0.02268114 28.47677       0
## 
##  Correlation: 
##         (Intr)
## log(BM) -0.297
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.91042516 -0.83939318 -0.05097508  0.63224435  3.05176705 
## 
## Residual standard error: 0.7847368 
## Degrees of freedom: 241 total; 239 residual
exp(-(5.050557/0.619536)+(1/0.619536)*log(5471405.789))
## [1] 21657214
data[1,8]<-21657214
dataV2<-subset(data, !is.na(V))

ggplot(dataV2, aes(log(BM), log(V), color=Suborder)) +
  geom_point(size=4) +
  xlab("ln body mass (g)") +
  ylab("ln egg volume (mm3)") +
  geom_abline(intercept=LambdaBM$coefficients[1], slope=LambdaBM$coefficients[2], colour="lightblue") +
  theme(panel.background = element_rect(fill="black")) +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15),
        legend.text=element_text(size=15),
        legend.title=element_text(size=17)) +
  scale_colour_brewer("Suborder", palette="Set1")

Eggshell thickness ~ Egg volume

  • Removing species with missing data
dataVT2<-subset(data2, !is.na(Thickness)&!is.na(V)); dataVT2<-dataVT2[-1,]
pruned.treeVT<-drop.tip(tree, setdiff(tree$tip.label, dataVT2$Species))
dataVT<-subset(data, !is.na(Thickness)&!is.na(V)); dataVT<-dataVT[-1,]
  • Building PGLS models and running model selection using AICc
BMVT<-gls(log(Thickness)~log(V), data=dataVT, correlation=corBrownian(phy=pruned.treeVT), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OUVT<-gls(log(Thickness)~log(V), data=dataVT, correlation=corMartins(1, phy=pruned.treeVT), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
LambdaVT<-gls(log(Thickness)~log(V), data=dataVT, correlation=corPagel(1, phy=pruned.treeVT), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Error in corFactor.corStruct(object): NA/NaN/Inf in foreign function call (arg 1)
EBVT<-gls(log(Thickness)~log(V), data=dataVT, correlation=corBlomberg(0.1, phy=pruned.treeVT, fixed=T), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OLSVT<-gls(log(Thickness)~log(V), data=dataVT, method="ML")

Cand.models = list()
Cand.models[[1]] = BMVT
Cand.models[[2]] = OUVT
Cand.models[[3]] = LambdaVT
## Error in eval(expr, envir, enclos): object 'LambdaVT' not found
Cand.models[[4]] = EBVT
Cand.models[[5]] = OLSVT

Modnames = paste(c("BMVT", "OUVT", "LambdaVT", "EBVT", "OLSVT"), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = T)
## Error in formatCands(cand.set): 
## Functions do not support mixture of model classes
summary(EBVT)
## Generalized least squares fit by maximum likelihood
##   Model: log(Thickness) ~ log(V) 
##   Data: dataVT 
##        AIC      BIC    logLik
##   147.5644 154.2229 -70.78221
## 
## Correlation Structure: corBlomberg
##  Formula: ~1 
##  Parameter estimate(s):
##   g 
## 0.1 
## 
## Coefficients:
##                 Value Std.Error  t-value p-value
## (Intercept) 2.3502382 0.4486431 5.238548       0
## log(V)      0.2840476 0.0555982 5.108931       0
## 
##  Correlation: 
##        (Intr)
## log(V) -0.966
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.6576055 -0.5886486 -0.1688905  0.5376755  2.7915997 
## 
## Residual standard error: 0.7289907 
## Degrees of freedom: 68 total; 66 residual
  • Testing for normality
lillie.test(residuals(EBVT))
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuals(EBVT)
## D = 0.071519, p-value = 0.5269
lillie.test(chol(solve(vcv(pruned.treeVT)))%*%residuals(EBVT))
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  chol(solve(vcv(pruned.treeVT))) %*% residuals(EBVT)
## D = 0.10736, p-value = 0.05007
plot(EBVT, resid(., type="n")~fitted(.), main="Normalized Residuals v Fitted Values",abline=c(0,0))

res<-resid(EBVT, type="n")
par(mar=c(5.1,4.1,4.1,2.1))
qqnorm(res)
qqline(res)

  • Removing outliers and building a new model without them
which(res<(-2)); which(res>2.5)
## Zootoca_vivipara_viviparous            Crotalus_viridis 
##                          45                          67
## Chamaeleo_senegalensis 
##                     19
newdataVT<-dataVT[-c(19,45),]
newdataVT2<-dataVT2[-c(19,45),]
newtreeVT<-drop.tip(pruned.treeVT, c("Zootoca_vivipara_viviparous","Chamaeleo_senegalensis"))

NewEBVT<-gls(log(Thickness)~log(V), data=newdataVT, correlation=corBlomberg(0.1, phy=newtreeVT, fixed=T), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
summary(NewEBVT)
## Generalized least squares fit by maximum likelihood
##   Model: log(Thickness) ~ log(V) 
##   Data: newdataVT 
##       AIC      BIC    logLik
##   126.345 132.9139 -60.17249
## 
## Correlation Structure: corBlomberg
##  Formula: ~1 
##  Parameter estimate(s):
##   g 
## 0.1 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 2.343769 0.3932038 5.960697       0
## log(V)      0.284315 0.0482883 5.887863       0
## 
##  Correlation: 
##        (Intr)
## log(V) -0.965
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.4180308 -0.6450608 -0.1851669  0.6097965  2.0464452 
## 
## Residual standard error: 0.6409687 
## Degrees of freedom: 66 total; 64 residual
plot(NewEBVT, resid(., type="n")~fitted(.), main="Normalized Residuals v Fitted Values",abline=c(0,0))

res<-resid(NewEBVT, type="n")
par(mar=c(5.1,4.1,4.1,2.1))
qqnorm(res)
qqline(res)

lillie.test(residuals(NewEBVT))
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuals(NewEBVT)
## D = 0.070589, p-value = 0.5719
lillie.test(chol(solve(vcv(newtreeVT)))%*%residuals(NewEBVT))
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  chol(solve(vcv(newtreeVT))) %*% residuals(NewEBVT)
## D = 0.088221, p-value = 0.2293
  • Plot of the regression
NewEBVT$coefficients
## (Intercept)      log(V) 
##    2.343769    0.284315
ggplot(newdataVT, aes(log(V), log(Thickness), color=Suborder)) +
  geom_point(size=4) +
  xlab("ln egg volume (mm3)") +
  ylab("ln eggshell thickness (µm)") +
  geom_abline(intercept=NewEBVT$coefficients[1], slope=NewEBVT$coefficients[2], colour="lightblue") +
  theme(panel.background = element_rect(fill="black")) +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15),
        legend.text=element_text(size=15),
        legend.title=element_text(size=17)) +
  scale_colour_brewer("Suborder", palette="Accent")

Phylogenetic path analysis

  • Building the models
ppadata<-dataVT
ppadata[,c(4,6:8)]<-log(ppadata[,c(4,6:8)])
M<-define_model_set(null=c(),
                    one=c(Thickness~SVL),
                    two=c(V~BM),
                    three=c(Thickness~V),
                    four=c(Thickness~V+SVL),
                    five=c(Thickness~V, V~BM),
                    six=c(Thickness~SVL, V~BM),
                    seven=c(Thickness~SVL+BM),
                    eight=c(V~BM, Thickness~BM),
                    nine=c(Thickness~V+BM),
                    ten=c(Thickness~SVL+BM+V),
                    eleven=c(Thickness~V+BM, V~BM),
                    twelve=c(Thickness~SVL+BM, V~BM),
                    thirteen=c(Thickness~BM),
                    .common=c(BM~SVL, V~SVL))
  • Plot all models and compile best model using CICc
plot_model_set(M)
## Warning: The curvature argument has been deprecated in favour of strength

result<-phylo_path(M, data=ppadata, tree=pruned.treeVT, model='lambda'); result
## A phylogenetic path analysis, on the variables:
##  Continuous:  Thickness SVL V BM 
##  Binary:       
## 
##  Evaluated for these models: null one two three four five six seven eight nine ten eleven twelve thirteen 
## 
##  Containing 31 phylogenetic regressions, of which 10 unique
s<-summary(result); s
plot(s)
## Warning: Position guide is perpendicular to the intended axis. Did you mean to
## specify a different guide `position`?

averagemodel<-average(result)
## Registered S3 methods overwritten by 'MuMIn':
##   method         from   
##   nobs.pgls      caper  
##   nobs.phylolm   phylolm
##   logLik.phylolm phylolm
plot(averagemodel)
## Warning: The curvature argument has been deprecated in favour of strength

Predictions of SVL and body mass – using phylogenetic eigenvectors maps (PEM)

  • Data and tree
fossiltreePEM<-read.nexus("Lepidosaurtree.trees.nex")
fossildataPEM<-read.table("Dataset1-lepidosaurs.txt",header=TRUE,stringsAsFactor=FALSE,)
fossildataV<-subset(fossildataPEM, !is.na(V))
fossildataV[,c(3:5,7:9)]<-log(fossildataV[,c(3:5,7:9)])

For SVL

  • Match the phylogeny and dataset
treePEM<-drop.tip(fossiltreePEM, setdiff(fossiltreePEM$tip.label, fossildataV$Species))
treedrop <- drop.tip(treePEM,fossildataV[is.na(fossildataV[,"SVL"]),"Species"])
grloc <- getGraphLocations(treePEM,fossildataV[is.na(fossildataV[,"SVL"]),"Species"])
sporder <- match(attr(grloc$x,"vlabel")[grloc$x$vertex$species],fossildataV[,"Species"])
  • Build the PEM
PEMfs <- list()
PEMfs[["V"]] <- PEM.fitSimple(y=fossildataV[sporder,"SVL"],
                    x=fossildataV[sporder,"V"],w=grloc$x,d="distance",sp="species",
                    lower=0,upper=1)
PEMfs[["none"]] <- PEM.fitSimple(y=fossildataV[sporder,"SVL"],
                    x=NULL,w=grloc$x,d="distance",sp="species",lower=0,upper=1)
for(m in c("V","none")) print(PEMfs[[m]]$optim$par)
## [1] 0.4981507
## [1] 0.1197493
rm(m)
  • Select the best model (with V as a co-predictor vs without co-predictors) based on AICc
PEMAIC <- list()
PEMAIC[["V"]] <- lmforwardsequentialAICc(y=fossildataV[sporder,"SVL"],
                                         x=fossildataV[sporder,"V",drop=FALSE],object=PEMfs[["V"]])
PEMAIC[["none"]] <- lmforwardsequentialAICc(y=fossildataV[sporder,"SVL"],object=PEMfs[["none"]])
for(m in c("V","none"))
  cat(m,summary(PEMAIC[[m]])$adj,PEMAIC[[m]]$AICc,"\n")
## V 0.9945104 -248.9384 
## none 0.9913577 -124.1351
rm(m)

summary(PEMAIC[["V"]])
## 
## Call:
## lm(formula = as.formula(paste(p1, p2, sep = "")), data = df1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.178636 -0.035423  0.004214  0.039172  0.153809 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.071635   0.066712  31.053  < 2e-16 ***
## V            0.393493   0.009214  42.705  < 2e-16 ***
## V_1         -9.100852   0.108193 -84.116  < 2e-16 ***
## V_2         -4.704946   0.126827 -37.097  < 2e-16 ***
## V_30         1.372985   0.083670  16.410  < 2e-16 ***
## V_6          0.986086   0.110409   8.931 3.32e-15 ***
## V_10         0.858355   0.118694   7.232 3.56e-11 ***
## V_8         -0.984514   0.083734 -11.758  < 2e-16 ***
## V_63        -0.760769   0.084320  -9.022 1.98e-15 ***
## V_7          0.862697   0.083710  10.306  < 2e-16 ***
## V_31         0.874262   0.084017  10.406  < 2e-16 ***
## V_117        0.819127   0.084435   9.701  < 2e-16 ***
## V_27        -1.240447   0.106374 -11.661  < 2e-16 ***
## V_77        -0.733867   0.083703  -8.768 8.30e-15 ***
## V_229       -0.819382   0.086023  -9.525  < 2e-16 ***
## V_156        0.507389   0.085388   5.942 2.39e-08 ***
## V_124       -0.408268   0.089180  -4.578 1.08e-05 ***
## V_5          0.388503   0.087559   4.437 1.91e-05 ***
## V_180       -0.572372   0.083671  -6.841 2.72e-10 ***
## V_14        -0.450315   0.085040  -5.295 4.87e-07 ***
## V_225       -0.546724   0.083745  -6.528 1.33e-09 ***
## V_48        -0.617529   0.083862  -7.364 1.77e-11 ***
## V_11         0.632745   0.084051   7.528 7.39e-12 ***
## V_108        0.675610   0.084757   7.971 6.76e-13 ***
## V_76         0.648948   0.084352   7.693 3.05e-12 ***
## V_44        -0.384803   0.085978  -4.476 1.64e-05 ***
## V_17         0.485665   0.083976   5.783 5.10e-08 ***
## V_80         0.535035   0.083675   6.394 2.61e-09 ***
## V_205        0.551436   0.083744   6.585 1.00e-09 ***
## V_23        -0.401879   0.084055  -4.781 4.61e-06 ***
## V_59        -0.320047   0.085352  -3.750 0.000265 ***
## V_66         0.544710   0.084217   6.468 1.81e-09 ***
## V_13        -0.499332   0.083867  -5.954 2.26e-08 ***
## V_29        -0.594629   0.085593  -6.947 1.57e-10 ***
## V_170        0.403913   0.083791   4.820 3.90e-06 ***
## V_109       -0.541708   0.084637  -6.400 2.54e-09 ***
## V_98        -0.558359   0.085053  -6.565 1.11e-09 ***
## V_179       -0.390290   0.083857  -4.654 7.86e-06 ***
## V_18        -0.529632   0.085127  -6.222 6.15e-09 ***
## V_206       -0.439153   0.083739  -5.244 6.13e-07 ***
## V_54        -0.453173   0.083883  -5.402 3.00e-07 ***
## V_125       -0.458915   0.084018  -5.462 2.28e-07 ***
## V_142        0.402447   0.083669   4.810 4.08e-06 ***
## V_235        0.429425   0.083804   5.124 1.05e-06 ***
## V_155        0.370097   0.083716   4.421 2.04e-05 ***
## V_34         0.351834   0.083814   4.198 4.94e-05 ***
## V_194        0.273009   0.085159   3.206 0.001692 ** 
## V_183       -0.376601   0.083670  -4.501 1.48e-05 ***
## V_169       -0.344993   0.083766  -4.119 6.71e-05 ***
## V_112        0.278725   0.084640   3.293 0.001274 ** 
## V_39        -0.293362   0.084304  -3.480 0.000682 ***
## V_55        -0.391209   0.083709  -4.673 7.26e-06 ***
## V_69         0.380311   0.083680   4.545 1.24e-05 ***
## V_129       -0.395590   0.083754  -4.723 5.89e-06 ***
## V_173       -0.381593   0.083692  -4.559 1.16e-05 ***
## V_152       -0.131500   0.088366  -1.488 0.139119    
## V_56         0.147425   0.087294   1.689 0.093629 .  
## V_62         0.424582   0.084164   5.045 1.48e-06 ***
## V_104        0.367887   0.083752   4.393 2.29e-05 ***
## V_72        -0.370855   0.083789  -4.426 2.00e-05 ***
## V_88        -0.321197   0.083667  -3.839 0.000191 ***
## V_178        0.381534   0.084033   4.540 1.26e-05 ***
## V_133       -0.435337   0.084932  -5.126 1.04e-06 ***
## V_58         0.361402   0.083871   4.309 3.19e-05 ***
## V_68         0.363490   0.084099   4.322 3.03e-05 ***
## V_172       -0.275813   0.083706  -3.295 0.001266 ** 
## V_176        0.264259   0.083676   3.158 0.001972 ** 
## V_151        0.262899   0.083679   3.142 0.002077 ** 
## V_101       -0.279064   0.083671  -3.335 0.001109 ** 
## V_130        0.354518   0.084439   4.199 4.93e-05 ***
## V_20         0.450383   0.087594   5.142 9.68e-07 ***
## V_16         0.364896   0.084960   4.295 3.38e-05 ***
## V_4          0.448939   0.088214   5.089 1.22e-06 ***
## V_134       -0.322932   0.084261  -3.833 0.000196 ***
## V_45        -0.265568   0.083671  -3.174 0.001874 ** 
## V_120        0.273162   0.083703   3.263 0.001404 ** 
## V_227       -0.273278   0.083759  -3.263 0.001408 ** 
## V_33         0.347907   0.085649   4.062 8.33e-05 ***
## V_148        0.267991   0.083752   3.200 0.001725 ** 
## V_118       -0.283059   0.083948  -3.372 0.000982 ***
## V_138        0.267968   0.083770   3.199 0.001731 ** 
## V_121       -0.253752   0.083678  -3.032 0.002925 ** 
## V_42         0.299744   0.084538   3.546 0.000544 ***
## V_97        -0.272391   0.083932  -3.245 0.001489 ** 
## V_116       -0.260448   0.083813  -3.107 0.002314 ** 
## V_73        -0.209929   0.083934  -2.501 0.013613 *  
## V_85         0.228299   0.083680   2.728 0.007242 ** 
## V_19         0.264901   0.084013   3.153 0.002003 ** 
## V_28        -0.307979   0.086077  -3.578 0.000486 ***
## V_26         0.224630   0.083682   2.684 0.008207 ** 
## V_21        -0.292578   0.085859  -3.408 0.000871 ***
## V_61        -0.242060   0.083830  -2.888 0.004544 ** 
## V_110       -0.242602   0.083878  -2.892 0.004479 ** 
## V_177        0.245034   0.084113   2.913 0.004208 ** 
## V_193       -0.215603   0.083688  -2.576 0.011095 *  
## V_99        -0.244967   0.085178  -2.876 0.004704 ** 
## V_95        -0.219343   0.084007  -2.611 0.010080 *  
## V_187        0.201047   0.083676   2.403 0.017677 *  
## V_115       -0.194921   0.083706  -2.329 0.021409 *  
## V_40        -0.204034   0.083748  -2.436 0.016183 *  
## V_35         0.231018   0.086910   2.658 0.008837 ** 
## V_38         0.223027   0.085936   2.595 0.010530 *  
## V_150       -0.203677   0.084038  -2.424 0.016732 *  
## V_222       -0.198031   0.083780  -2.364 0.019564 *  
## V_168       -0.188724   0.083980  -2.247 0.026296 *  
## V_141        0.183297   0.083777   2.188 0.030449 *  
## V_12         0.180644   0.083706   2.158 0.032745 *  
## V_91        -0.174506   0.084069  -2.076 0.039873 *  
## V_201       -0.168706   0.083691  -2.016 0.045863 *  
## V_67        -0.162550   0.083667  -1.943 0.054182 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08367 on 131 degrees of freedom
## Multiple R-squared:  0.997,  Adjusted R-squared:  0.9945 
## F-statistic: 399.9 on 109 and 131 DF,  p-value: < 2.2e-16

The model with V as a co-predictor has the best fit.

  • Predicting missing values for SVL
tf<-list(V=identity, none=identity)
m <- "V" ; atr <- tf[[m]](fossildataV[is.na(fossildataV[,"SVL"]),m,drop=FALSE])
resultsPEM<-predict(object=PEMfs[[m]],targets=grloc,lmobject=PEMAIC[[m]],newdata=atr,interval="confidence")

# Predicted SVL, with upper and lower limits of the confidence interval for the prediction:
exp(resultsPEM$values); exp(resultsPEM$upper); exp(resultsPEM$lower)
##                              [,1]
## Antarcticoolithus_bradyi 6665.867
##                              [,1]
## Antarcticoolithus_bradyi 7888.606
##                              [,1]
## Antarcticoolithus_bradyi 5632.654

The code to predict missing values for body mass is identical to the one for SVL; simply replace ‘SVL’ with ‘BM’.

  • Plotting the regression for Egg volume ~ SVL with the value of SVL for Antarcticoolithus estimated from PEM
data[1,7]<-6665.867
dataV2<-subset(data, !is.na(V))

ggplot(dataV2, aes(log(SVL), log(V), color=Suborder)) +
  geom_point(size=4) +
  xlab("ln snout-vent length (mm)") +
  ylab("ln egg volume (mm3)") +
  geom_abline(intercept=Lambda$coefficients[1], slope=Lambda$coefficients[2], colour="lightblue") +
  theme(panel.background = element_rect(fill="black")) +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15),
        legend.text=element_text(size=15),
        legend.title=element_text(size=17)) +
  scale_colour_brewer("Suborder", palette="Accent")

  • Plotting the regression for Egg volume ~ Body mass with the value of body mass for Antarcticoolithus estimated from PEM
data[1,8]<-456096.7
dataV2<-subset(data, !is.na(V))

ggplot(dataV2, aes(log(BM), log(V), color=Suborder)) +
  geom_point(size=4) +
  xlab("ln body mass (g)") +
  ylab("ln egg volume (mm3)") +
  geom_abline(intercept=LambdaBM$coefficients[1], slope=LambdaBM$coefficients[2], colour="lightblue") +
  theme(panel.background = element_rect(fill="black")) +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15),
        legend.text=element_text(size=15),
        legend.title=element_text(size=17)) +
  scale_colour_brewer("Suborder", palette="Set1")

Phylogenetic ANOVA

  • Subsetting the data and tree, and compiling relative eggshell thickness (RT) with RT = (eggshell thickness / egg length)
RT<-data$Thickness/data$L; data<-cbind(data, RT)
dataA<-subset(data, !is.na(RT)&!is.na(OV))
dataA$RT<-log(dataA$RT)
anovatree<-drop.tip(tree, setdiff(tree$tip.label, rownames(dataA)))

ANOVA with 2 groups (oviparous and viviparous)

  • Subsetting the dataset and performing the ANOVA
datANOVA<-as.data.frame(cbind(dataA$RT, dataA$OV), header=T)
row.names(datANOVA)<-row.names(dataA); colnames(datANOVA)<-c("RT", "OV")
datANOVA<-ReorderData(anovatree, datANOVA, taxa.names="row names")

phylANOVA(anovatree, datANOVA$OV, datANOVA$RT, p.adj="BH")
## Warning: no labels for x. Assuming order of tree$tip.label.
## 
## Warning: no labels for y. Assuming order of tree$tip.label.
## ANOVA table: Phylogenetic ANOVA
## 
## Response: y
##             Sum Sq  Mean Sq  F value Pr(>F)
## x         2.884385 2.884385 6.074072  0.011
## Residual 31.816186 0.474868                
## 
## P-value based on simulation.
## ---------
## 
## Pairwise posthoc test using method = "BH"
## 
## Pairwise t-values:
##           1        2
## 1  0.000000 2.464563
## 2 -2.464563 0.000000
## 
## Pairwise corrected P-values:
##       1     2
## 1 1.000 0.011
## 2 0.011 1.000
## ---------
  • Boxplot of both groups
ggplot(dataA, aes(x=OV, y=RT, fill=OV)) +
  geom_boxplot() +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15),
        legend.text=element_text(size=15),
        legend.title=element_text(size=17)) +
  scale_fill_hue()

log(data$RT)[1] # Value for Antarcticoolithus, intermediate between the groups
## [1] 0.9120565
  • Extract percentiles for each boxplot
p2<-ggplot(dataA, aes(x=OV, y=RT, fill=OV)) + geom_boxplot()
ggplot_build(p2)$data
## [[1]]
##      fill       ymin      lower      middle     upper     ymax outliers
## 1 #F8766D -0.1251631  0.9808293  1.35720692 1.7429693 2.784423 3.442019
## 2 #00BFC4 -0.3462762 -0.3308263 -0.02868824 0.8270896 2.503459         
##   notchupper notchlower x flipped_aes PANEL group ymin_final ymax_final  xmin
## 1  1.5065672  1.2078466 1       FALSE     1     1 -0.1251631   3.442019 0.625
## 2  0.8860653 -0.9434418 2       FALSE     1     2 -0.3462762   2.503459 1.625
##    xmax xid newx new_width weight colour size alpha shape linetype
## 1 1.375   1    1      0.75      1 grey20  0.5    NA    19    solid
## 2 2.375   2    2      0.75      1 grey20  0.5    NA    19    solid

ANOVA with 3 groups (oviparous and laid, oviparous and dissected from the oviduct, and viviparous)

  • Subsetting the dataset and performing the ANOVA
datANOVA2<-as.data.frame(cbind(dataA$RT, dataA$Strategy), header=T)
row.names(datANOVA2)<-row.names(dataA); colnames(datANOVA2)<-c("RT", "Strategy")
datANOVA2<-ReorderData(anovatree, datANOVA2, taxa.names="row names")

phylANOVA(anovatree, datANOVA2$Strategy, datANOVA2$RT, p.adj="BH")
## Warning: no labels for x. Assuming order of tree$tip.label.
## 
## Warning: no labels for y. Assuming order of tree$tip.label.
## ANOVA table: Phylogenetic ANOVA
## 
## Response: y
##            Sum Sq  Mean Sq  F value Pr(>F)
## x         7.86154 3.930770 9.666176  0.005
## Residual 26.83903 0.406652                
## 
## P-value based on simulation.
## ---------
## 
## Pairwise posthoc test using method = "BH"
## 
## Pairwise t-values:
##           1         2        3
## 1  0.000000 -1.534310 3.960211
## 2  1.534310  0.000000 4.382854
## 3 -3.960211 -4.382854 0.000000
## 
## Pairwise corrected P-values:
##       1     2     3
## 1 1.000 0.311 0.003
## 2 0.311 1.000 0.003
## 3 0.003 0.003 1.000
## ---------
  • Boxplot of all groups
ggplot(dataA, aes(x=Strategy, y=RT, fill=Strategy)) +
  geom_boxplot() +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15),
        legend.text=element_text(size=15),
        legend.title=element_text(size=17)) +
  scale_fill_hue()

  • Extract percentiles for each boxplot
p2<-ggplot(dataA, aes(x=Strategy, y=RT, fill=Strategy)) + geom_boxplot()
ggplot_build(p2)$data
## [[1]]
##      fill       ymin      lower     middle       upper      ymax outliers
## 1 #F8766D -0.1251631  0.9794977  1.3139737  1.74296931 2.7844232         
## 2 #00BA38  0.3137979  1.3231783  1.6757606  2.06667501 2.5034588 3.442019
## 3 #619CFF -0.3462762 -0.3359763 -0.3256763 -0.02868824 0.2682998         
##    notchupper notchlower x flipped_aes PANEL group ymin_final ymax_final  xmin
## 1  1.47662916  1.1513182 1       FALSE     1     1 -0.1251631  2.7844232 0.625
## 2  2.02995340  1.3215677 2       FALSE     1     2  0.3137979  3.4420194 1.625
## 3 -0.04536403 -0.6059886 3       FALSE     1     3 -0.3462762  0.2682998 2.625
##    xmax xid newx new_width weight colour size alpha shape linetype
## 1 1.375   1    1      0.75      1 grey20  0.5    NA    19    solid
## 2 2.375   2    2      0.75      1 grey20  0.5    NA    19    solid
## 3 3.375   3    3      0.75      1 grey20  0.5    NA    19    solid

2 – Amniote dataset

This part requires the use of dataset 2 (“Dataset2-amniotes.txt”; see also Supplementary Table 3 in Legendre et al., 2020), and amniote phylogenetic tree (“Amniotetree.nex”).

  • Loading and matching the data and tree
dataS<-read.table("Dataset2-amniotes.txt", header=T)
treeS<-read.nexus("Amniotetree.nex")
treeS<-drop.tip(treeS, setdiff(treeS$tip.label, dataS$Taxon))
dataS<-read.table("Dataset2-amniotes.txt", header=T, row.names="Taxon")
Ws<-diag(vcv.phylo(treeS))

PGLS – Eggshell thickness ~ Egg mass

For hard-shelled eggs

  • Subset data and tree
dataS<-read.table("Dataset2-amniotes.txt", header=T)
datahard<-dataS[-c(50,86:88,93:95,115,118,120:145,147,149),]
treehard<-drop.tip(treeS, setdiff(treeS$tip.label, datahard$Taxon))
dataS<-read.table("Dataset2-amniotes.txt", header=T, row.names="Taxon")
datahard<-dataS[-c(50,86:88,93:95,115,118,120:145,147,149),]
Wh<-diag(vcv.phylo(treehard))
  • Best fit for alpha parameter in OU model
alpha <- seq(0, 1, 0.1)
fit <- list()
form <- log(Eggshell_thickness)~log(Egg_mass)
for (i in seq_along(alpha)) {
  cor <- corMartins(alpha[i], phy = treehard, fixed = T)
  fit[[i]] <- gls(form, correlation = cor, data = datahard, na.action=na.exclude, weights=varFixed(~Wh), method = "ML")
}
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
sapply(fit, logLik)
##  [1]       Inf -37.10514 -37.00214 -37.38638 -37.59015 -37.68151 -37.71941
##  [8] -37.73459 -37.74057 -37.74291 -37.74382
  • Best fit for g parameter in EB (Early Burst) model
g <- seq(0.1, 1, 0.1)
fit <- list()
form <- log(Eggshell_thickness)~log(Egg_mass)
for (i in seq_along(g)) {
  cor <- corBlomberg(g[i], phy = treehard, fixed = T)
  fit[[i]] <- gls(form, correlation = cor, data = datahard, na.action=na.exclude, weights=varFixed(~Wh), method = "ML")
}
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
sapply(fit, logLik)
##  [1] -51.97396 -62.61195 -68.20028 -71.62168 -73.95198 -75.65088 -76.94873
##  [8] -77.97487 -78.80802 -79.49898
  • Building PGLS models and running model selection using AICc
BM<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datahard, correlation=corBrownian(phy=treehard), weights=varFixed(~Wh), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OU<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datahard, correlation=corMartins(0.1, phy=treehard, fixed=T), weights=varFixed(~Wh), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
Lambda<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datahard, correlation=corPagel(1, phy=treehard), weights=varFixed(~Wh), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
EB<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datahard, correlation=corBlomberg(0.1, phy=treehard, fixed=T), weights=varFixed(~Wh), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OLS<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datahard, method="ML")

Cand.models = list()
Cand.models[[1]] = BM
Cand.models[[2]] = OU
Cand.models[[3]] = Lambda
Cand.models[[4]] = EB
Cand.models[[5]] = OLS

Modnames = paste(c("BM", "OU", "Lambda", "EB", "OLS"), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = T)
summary(Lambda)
## Generalized least squares fit by maximum likelihood
##   Model: log(Eggshell_thickness) ~ log(Egg_mass) 
##   Data: datahard 
##        AIC      BIC    logLik
##   67.05542 77.92942 -29.52771
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.3659904 
## Variance function:
##  Structure: fixed weights
##  Formula: ~Wh 
## 
## Coefficients:
##                  Value  Std.Error  t-value p-value
## (Intercept)   4.386080 0.12772437 34.34020       0
## log(Egg_mass) 0.426746 0.01576025 27.07737       0
## 
##  Correlation: 
##               (Intr)
## log(Egg_mass) -0.443
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.78157887 -0.99474516 -0.57993810  0.03130181  2.28401650 
## 
## Residual standard error: 0.02035964 
## Degrees of freedom: 112 total; 110 residual
  • Estimating the pseudo R-squared using caper
dataS<-read.table("Dataset2-amniotes.txt", header=T)
datahard<-dataS[-c(50,86:88,93:95,115,118,120:145,147,149),]
datacomp<-comparative.data(phy=treehard, data=datahard, names.col="Taxon")
pgls<-pgls(log(Eggshell_thickness)~log(Egg_mass), data=datacomp, lambda="ML")
summary(pgls)
## 
## Call:
## pgls(formula = log(Eggshell_thickness) ~ log(Egg_mass), data = datacomp, 
##     lambda = "ML")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.06443 -0.01539 -0.00283  0.01159  0.05560 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.328
##    lower bound : 0.000, p = 0.0012278
##    upper bound : 1.000, p = < 2.22e-16
##    95.0% CI   : (0.098, 0.597)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   4.334504   0.128404  33.757 < 2.2e-16 ***
## log(Egg_mass) 0.436835   0.017409  25.092 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02103 on 110 degrees of freedom
## Multiple R-squared: 0.8513,  Adjusted R-squared: 0.8499 
## F-statistic: 629.6 on 1 and 110 DF,  p-value: < 2.2e-16

For soft-shelled eggs

  • Subset data and tree
datasoft<-dataS[c(50,86:88,93,94,115,118,120:145,147,149),]
treesoft<-drop.tip(treeS, setdiff(treeS$tip.label, datasoft$Taxon))
dataS<-read.table("Dataset2-amniotes.txt", header=T, row.names="Taxon")
datasoft<-dataS[c(50,86:88,93,94,115,118,120:145,147,149),]
Ws<-diag(vcv.phylo(treesoft))
  • Best fit for alpha parameter in OU model
alpha <- seq(0, 1, 0.1)
fit <- list()
form <- log(Eggshell_thickness)~log(Egg_mass)
for (i in seq_along(alpha)) {
  cor <- corMartins(alpha[i], phy = treesoft, fixed = T)
  fit[[i]] <- gls(form, correlation = cor, data = datasoft, na.action=na.exclude, weights=varFixed(~Ws), method = "ML")
}
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
sapply(fit, logLik)
##  [1]       Inf -31.45442 -31.58754 -31.70751 -31.76148 -31.78240 -31.79013
##  [8] -31.79293 -31.79394 -31.79430 -31.79443
  • Best fit for g parameter in EB model
g <- seq(0.1, 1, 0.1)
fit <- list()
form <- log(Eggshell_thickness)~log(Egg_mass)
for (i in seq_along(g)) {
  cor <- corBlomberg(g[i], phy = treesoft, fixed = T)
  fit[[i]] <- gls(form, correlation = cor, data = datasoft, na.action=na.exclude, weights=varFixed(~Ws), method = "ML")
}
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.

## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
sapply(fit, logLik)
##  [1] -31.02830 -33.27681 -35.31325 -36.86832 -38.06022 -38.99694 -39.75253
##  [8] -40.37614 -40.90094 -41.34988
  • Building PGLS models and running model selection using AICc
BM<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datasoft, correlation=corBrownian(phy=treesoft), weights=varFixed(~Ws), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OU<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datasoft, correlation=corMartins(0.2, phy=treesoft), weights=varFixed(~Ws), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Error in recalc.corStruct(object[[i]], conLin): NA/NaN/Inf in foreign function call (arg 4)
Lambda<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datasoft, correlation=corPagel(1, phy=treesoft), weights=varFixed(~Ws), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
EB<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datasoft, correlation=corBlomberg(0.1, phy=treesoft, fixed=T), weights=varFixed(~Ws), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OLS<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datasoft, method="ML")

Cand.models = list()
Cand.models[[1]] = BM
Cand.models[[2]] = OU
Cand.models[[3]] = Lambda
Cand.models[[4]] = EB
Cand.models[[5]] = OLS

Modnames = paste(c("BM", "OU", "Lambda", "EB", "OLS"), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = T)
summary(EB)
## Generalized least squares fit by maximum likelihood
##   Model: log(Eggshell_thickness) ~ log(Egg_mass) 
##   Data: datasoft 
##        AIC      BIC   logLik
##   68.05659 72.80715 -31.0283
## 
## Correlation Structure: corBlomberg
##  Formula: ~1 
##  Parameter estimate(s):
##   g 
## 0.1 
## Variance function:
##  Structure: fixed weights
##  Formula: ~Ws 
## 
## Coefficients:
##                   Value  Std.Error   t-value p-value
## (Intercept)   2.0167850 0.16983372 11.875056       0
## log(Egg_mass) 0.4009087 0.04963276  8.077501       0
## 
##  Correlation: 
##               (Intr)
## log(Egg_mass) -0.643
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.0217026 -0.6275853 -0.1318103  0.6864130  1.6585487 
## 
## Residual standard error: 0.03771851 
## Degrees of freedom: 36 total; 34 residual
  • Estimating the pseudo R-squared by compiling an ANOVA with the EB model and a null model
# R-squared
EB2<-gls(log(Eggshell_thickness)~1, data=datasoft, correlation=corBlomberg(0.1, phy=treesoft, fixed=T), weights=varFixed(~Ws), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
1 - (EB$sigma/EB2$sigma)^2
## [1] 0.657417
# P-value
anova(EB2, EB)

Plot with both regression lines (hard-shelled and soft-shelled)

dataS<-read.table("Dataset2-amniotes.txt", header=T)
ggplot(dataS, aes(log(Egg_mass), log(Eggshell_thickness), colour=Clade)) +
  geom_point(size=4) +
  xlab("ln egg mass (g)") +
  ylab("ln calcareous layer thickness (µm)") +
  geom_abline(intercept=4.321066, slope=0.439690, colour="turquoise", size=1.3) +
  geom_abline(intercept=2.0474668, slope=0.4036309, colour="red", size=1.3) +
  theme(panel.background = element_rect(fill="black")) +
  theme(axis.text=element_text(size=15),
        axis.title=element_text(size=15),
        legend.text=element_text(size=15),
        legend.title=element_text(size=17)) +
  scale_colour_brewer("Group", palette="Set1")

Ancestral state reconstruction

  • For ratio (Calcareous layer thickness / Egg mass), using fastAnc
dataS<-read.table("Dataset2-amniotes.txt", header=T); dataS<-dataS[-95,]
treeplot<-drop.tip(treeS, setdiff(treeS$tip.label, dataS$Taxon))
dataS<-read.table("Dataset2-amniotes.txt", header=T, row.names="Taxon"); dataS<-dataS[-95,]
dataplot=as.matrix(dataS[,3]); names(dataplot)<-rownames(dataS); colnames(dataplot)<-"CL_thickness"
dataplot<-log10(dataplot)

fit<-fastAnc(treeplot,dataplot,vars=T,CI=T) # Ancestral states for each node
obj<-contMap(treeplot,dataplot,plot=F)

plot(setMap(obj, colors=rev(brewer.pal(10, "Spectral"))), fsize=0.4,lwd=4)

  • Same as above, but with an ultrametric tree, as shown in Fig. 3 of Legendre et al. (2020)
treeplot<-force.ultrametric(treeplot)
fit<-fastAnc(treeplot,dataplot,vars=T,CI=T) # Ancestral states for each node
obj<-contMap(treeplot,dataplot,plot=F)

plot(setMap(obj, colors=rev(brewer.pal(10, "Spectral"))), fsize=0.4,lwd=4)

  • For presence/absence of the prismatic layer, using make.simmap
dataS<-read.table("Dataset2-amniotes.txt", header=T)
dataS<-dataS[-95,]
treeS<-read.nexus("Amniotetree.nex")
treeS<-drop.tip(treeS, setdiff(treeS$tip.label, dataS$Taxon))
prismatic<-dataS[,5]; names(prismatic)<-dataS$Taxon
cols<-setNames(c("royalblue","red3"),levels(prismatic))

mtree<-make.simmap(treeS, prismatic, model="ER")
## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##                Absent       Present
## Absent  -0.0003949281  0.0003949281
## Present  0.0003949281 -0.0003949281
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##  Absent Present 
##     0.5     0.5
## Done.
plot(mtree,cols,fsize=0.5,ftype="i")
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
                  y=0.8*par()$usr[3],fsize=0.8)

# with 1000 simulations of stochastic character maps from the data:
mtrees<-make.simmap(treeS, prismatic, model="ER", nsim=1000)
## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##                Absent       Present
## Absent  -0.0003949281  0.0003949281
## Present  0.0003949281 -0.0003949281
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##  Absent Present 
##     0.5     0.5
## Done.
par(mfrow=c(10,10))
null<-sapply(mtrees,plotSimmap,colors=cols,lwd=1,ftype="off")

pd<-summary(mtrees); pd
## 1000 trees with a mapped discrete character with states:
##  Absent, Present 
## 
## trees have 3.735 changes between states on average
## 
## changes are of the following types:
##      Absent,Present Present,Absent
## x->y          0.981          2.754
## 
## mean total time spent in each state is:
##            Absent      Present   total
## raw  2705.1756568 6756.0840379 9461.26
## prop    0.2859213    0.7140787    1.00
par(mfrow=c(1,1))
plot(pd,fsize=0.6,ftype="i",colors=cols,ylim=c(-2,Ntip(treeS)))
add.simmap.legend(colors=cols[2:1],prompt=FALSE,x=0,y=-2,vertical=FALSE)

# density map of the above
objprism<-densityMap(mtrees,states=levels(prismatic)[2:1],plot=FALSE)
## sorry - this might take a while; please be patient
n<-length(objprism$cols)
objprism$cols[1:n]<-colorRampPalette(c("royalblue","red3"))(n)
plot(objprism,size=c(0.6,1),fsize=0.6,lwd=4)

Boxplot for ratio (Calcareous layer thickness / Egg mass)

dataS<-read.table("Dataset2-amniotes.txt", header=T)
RTM<-dataS$Eggshell_thickness/dataS$Egg_mass
dataS<-cbind(dataS, RTM)
dataSboxplot<-dataS %>% filter(!Clade=="Test")

ggplot(dataSboxplot, aes(x=Clade, y=log(RTM), fill=Clade)) + 
  geom_boxplot() +
  theme(panel.background = element_rect(fill="white")) +
  theme(axis.text=element_text(size=10, angle = 90),
        axis.title=element_text(size=15),
        legend.text=element_text(size=15),
        legend.title=element_text(size=17)) +
  scale_fill_brewer(palette="Set1")

  • Extract percentiles for each boxplot
p<-ggplot(dataSboxplot, aes(x=Clade, y=log(RTM), fill=Clade)) + 
  geom_boxplot()
ggplot_build(p)$data
## [[1]]
##      fill       ymin      lower     middle      upper       ymax
## 1 #F8766D -4.1737717 -4.1737717 -4.1737717 -4.1737717 -4.1737717
## 2 #D39200  0.4296871  1.8265428  2.4120341  2.9394662  4.1207778
## 3 #93AA00  0.2636401  1.6216397  2.5162801  2.5794171  2.9038769
## 4 #00BA38  1.4714471  1.4714471  1.4714471  1.4714471  1.4714471
## 5 #00C19F  1.1282254  1.5449861  1.7324607  1.9034789  2.2137481
## 6 #00B9E3  4.4770380  4.4770380  4.4770380  4.4770380  4.4770380
## 7 #619CFF -2.1102132 -0.4793495 -0.1364957  0.9978001  1.9902104
## 8 #DB72FB  0.2876821  0.3729137  0.4581454  0.5433770  0.6286087
## 9 #FF61C3 -1.4307461  0.2392498  1.7701329  3.4510661  5.4647026
##                 outliers notchupper notchlower x flipped_aes         y PANEL
## 1                        -4.1737717 -4.1737717 1       FALSE -4.173772     1
## 2 -0.9350988, -0.6852781  2.6607121  2.1633561 2       FALSE        NA     1
## 3                         3.0207095  2.0118506 3       FALSE        NA     1
## 4                         1.4714471  1.4714471 4       FALSE  1.471447     1
## 5                         1.9327199  1.5322015 5       FALSE        NA     1
## 6                         4.4770380  4.4770380 6       FALSE  4.477038     1
## 7                         0.2372268 -0.5102183 7       FALSE        NA     1
## 8                         0.6485919  0.2676989 8       FALSE        NA     1
## 9                         2.5933530  0.9469128 9       FALSE        NA     1
##   group ymin_final ymax_final  xmin  xmax xid newx new_width weight colour size
## 1     1 -4.1737717 -4.1737717 0.625 1.375   1    1      0.75      1 grey20  0.5
## 2     2 -0.9350988  4.1207778 1.625 2.375   2    2      0.75      1 grey20  0.5
## 3     3  0.2636401  2.9038769 2.625 3.375   3    3      0.75      1 grey20  0.5
## 4     4  1.4714471  1.4714471 3.625 4.375   4    4      0.75      1 grey20  0.5
## 5     5  1.1282254  2.2137481 4.625 5.375   5    5      0.75      1 grey20  0.5
## 6     6  4.4770380  4.4770380 5.625 6.375   6    6      0.75      1 grey20  0.5
## 7     7 -2.1102132  1.9902104 6.625 7.375   7    7      0.75      1 grey20  0.5
## 8     8  0.2876821  0.6286087 7.625 8.375   8    8      0.75      1 grey20  0.5
## 9     9 -1.4307461  5.4647026 8.625 9.375   9    9      0.75      1 grey20  0.5
##   alpha shape linetype
## 1    NA    19    solid
## 2    NA    19    solid
## 3    NA    19    solid
## 4    NA    19    solid
## 5    NA    19    solid
## 6    NA    19    solid
## 7    NA    19    solid
## 8    NA    19    solid
## 9    NA    19    solid